There are a variety of helper functions. Check out some of the helper files: ## spatial correlation
?spat_cor_mat
example(spat_cor_mat)##
## spt_c_> #a correlation matrix with 10 sites and sd of 4
## spt_c_> round(spat_cor_mat(n_loc=10,spat_sd=4,sumto1=FALSE),3)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1.000 0.898 0.871 0.838 0.812 0.802 0.812 0.838 0.871 0.898
## [2,] 0.898 1.000 0.898 0.871 0.838 0.812 0.802 0.812 0.838 0.871
## [3,] 0.871 0.898 1.000 0.898 0.871 0.838 0.812 0.802 0.812 0.838
## [4,] 0.838 0.871 0.898 1.000 0.898 0.871 0.838 0.812 0.802 0.812
## [5,] 0.812 0.838 0.871 0.898 1.000 0.898 0.871 0.838 0.812 0.802
## [6,] 0.802 0.812 0.838 0.871 0.898 1.000 0.898 0.871 0.838 0.812
## [7,] 0.812 0.802 0.812 0.838 0.871 0.898 1.000 0.898 0.871 0.838
## [8,] 0.838 0.812 0.802 0.812 0.838 0.871 0.898 1.000 0.898 0.871
## [9,] 0.871 0.838 0.812 0.802 0.812 0.838 0.871 0.898 1.000 0.898
## [10,] 0.898 0.871 0.838 0.812 0.802 0.812 0.838 0.871 0.898 1.000
##
## spt_c_> #a stray th 10 sites and sd of 4
## spt_c_> round(spat_cor_mat(n_loc=10,spat_sd=4,sumto1=TRUE),3)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 0.116 0.104 0.101 0.097 0.094 0.093 0.094 0.097 0.101 0.104
## [2,] 0.104 0.116 0.104 0.101 0.097 0.094 0.093 0.094 0.097 0.101
## [3,] 0.101 0.104 0.116 0.104 0.101 0.097 0.094 0.093 0.094 0.097
## [4,] 0.097 0.101 0.104 0.116 0.104 0.101 0.097 0.094 0.093 0.094
## [5,] 0.094 0.097 0.101 0.104 0.116 0.104 0.101 0.097 0.094 0.093
## [6,] 0.093 0.094 0.097 0.101 0.104 0.116 0.104 0.101 0.097 0.094
## [7,] 0.094 0.093 0.094 0.097 0.101 0.104 0.116 0.104 0.101 0.097
## [8,] 0.097 0.094 0.093 0.094 0.097 0.101 0.104 0.116 0.104 0.101
## [9,] 0.101 0.097 0.094 0.093 0.094 0.097 0.101 0.104 0.116 0.104
## [10,] 0.104 0.101 0.097 0.094 0.093 0.094 0.097 0.101 0.104 0.116
##
## spt_c_> # illustrate different parameter values
## spt_c_> library(RColorBrewer)
##
## spt_c_> a <- data.frame(n_loc= rep(10,3),spat_scale=c(1,4,10),spat_sd=c(1,4,10))
##
## spt_c_> b <- array(cbind(with(a,mapply(spat_cor_mat,spat_scale=spat_scale,n_loc=n_loc)),
## spt_c_+ with(a,mapply(spat_cor_mat,spat_sd=spat_sd,n_loc=n_loc))),dim= c(10,10,6))
##
## spt_c_> df1 <- melt(b,varnames = c("site1","site2","parameter"),value.name="correlation")
##
## spt_c_> df1$site1 <- factor(df1$site1)
##
## spt_c_> df1$parameter <- factor(df1$parameter, labels=
## spt_c_+ c(paste(rep(c("Cauchy, scale=","Gaussian, sd="),each= 3),rep(c(1,4,10),2))))
##
## spt_c_> ggplot(aes(site2,correlation),data =df1)+
## spt_c_+ geom_line(aes(colour=site1))+
## spt_c_+ labs(colour="site")+
## spt_c_+ scale_x_continuous(breaks= c(1:10))+
## spt_c_+ xlab("site")+
## spt_c_+ facet_wrap(~parameter)+theme_bw()
## spatiotemporal correlation
?spat_temp_ts
example(spat_temp_ts)##
## spt_t_> ### generate the time series
## spt_t_> errors <- spat_temp_ts(n_iter = 100, n_loc = 5,
## spt_t_+ site_sd = 0.5, spat_sd = 2,
## spt_t_+ phi = 0.5)
##
## spt_t_> # examine correlation matrix of time series
## spt_t_> errors$cor_mat
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0000000 0.8260514 0.5969475 0.6017253 0.7671826
## [2,] 0.8260514 1.0000000 0.8163391 0.6186857 0.5632422
## [3,] 0.5969475 0.8163391 1.0000000 0.7765027 0.5311337
## [4,] 0.6017253 0.6186857 0.7765027 1.0000000 0.7644076
## [5,] 0.7671826 0.5632422 0.5311337 0.7644076 1.0000000
##
## spt_t_> # examine pacf of each time series
## spt_t_> errors$pacf
## [1] 0.4938427 0.5033885 0.5299877 0.5191227 0.4103612
##
## spt_t_> # examine sd of each time series
## spt_t_> apply(errors$ts,2,sd)
## [1] 0.5630387 0.5339049 0.5269409 0.5730403 0.4730050
##
## spt_t_> # plot the time series
## spt_t_> matplot(errors$ts,type= "l")
?ran_stray_prob
example(ran_stray_prob)##
## rn_st_> stray_mat <- spat_cor_mat(n_loc=10,spat_sd=2, sumto1=TRUE)
##
## rn_st_> stray_mat
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.13732540 0.11902069 0.10503539 0.08999813 0.07942307 0.07572006
## [2,] 0.11902069 0.13732540 0.11902069 0.10503539 0.08999813 0.07942307
## [3,] 0.10503539 0.11902069 0.13732540 0.11902069 0.10503539 0.08999813
## [4,] 0.08999813 0.10503539 0.11902069 0.13732540 0.11902069 0.10503539
## [5,] 0.07942307 0.08999813 0.10503539 0.11902069 0.13732540 0.11902069
## [6,] 0.07572006 0.07942307 0.08999813 0.10503539 0.11902069 0.13732540
## [7,] 0.07942307 0.07572006 0.07942307 0.08999813 0.10503539 0.11902069
## [8,] 0.08999813 0.07942307 0.07572006 0.07942307 0.08999813 0.10503539
## [9,] 0.10503539 0.08999813 0.07942307 0.07572006 0.07942307 0.08999813
## [10,] 0.11902069 0.10503539 0.08999813 0.07942307 0.07572006 0.07942307
## [,7] [,8] [,9] [,10]
## [1,] 0.07942307 0.08999813 0.10503539 0.11902069
## [2,] 0.07572006 0.07942307 0.08999813 0.10503539
## [3,] 0.07942307 0.07572006 0.07942307 0.08999813
## [4,] 0.08999813 0.07942307 0.07572006 0.07942307
## [5,] 0.10503539 0.08999813 0.07942307 0.07572006
## [6,] 0.11902069 0.10503539 0.08999813 0.07942307
## [7,] 0.13732540 0.11902069 0.10503539 0.08999813
## [8,] 0.11902069 0.13732540 0.11902069 0.10503539
## [9,] 0.10503539 0.11902069 0.13732540 0.11902069
## [10,] 0.08999813 0.10503539 0.11902069 0.13732540
##
## rn_st_> rand_strays1 <- ran_stray_prob(stray_mat,n_iter=3,scale=100)
##
## rn_st_> rand_strays2 <- ran_stray_prob(stray_mat,n_iter=3,scale=1000)
##
## rn_st_> rand_strays3 <- ran_stray_prob(stray_mat,n_iter=3,scale=10^5)
##
## rn_st_> df1 <- melt(rand_strays1,varnames = c("site1","site2","replicate"),value.name="stray_probability")
##
## rn_st_> df2 <- melt(rand_strays2,varnames = c("site1","site2","replicate"),value.name="stray_probability")
##
## rn_st_> df3 <- melt(rand_strays3,varnames = c("site1","site2","replicate"),value.name="stray_probability")
##
## rn_st_> df4 <- rbind(df1,df2,df3)
##
## rn_st_> df4$precision <- factor(rep(c(1,2,3),each=nrow(df4)/3), labels= c("Precision = 10^2", "Precision =10^3", "Precision =10^5"))
##
## rn_st_> df4$site1 <- factor(df4$site1)
##
## rn_st_> ggplot(aes(site2,stray_probability),data =df4)+
## rn_st_+ geom_line(aes(colour=site1))+
## rn_st_+ labs(colour="source site")+ylab("straying probability (random Dirichlet samples")+
## rn_st_+ scale_x_continuous(breaks= c(1:10))+
## rn_st_+ xlab("destination site")+
## rn_st_+ facet_grid(precision~replicate)+theme_bw()
Run a simulation where there are 10 stages (1 egg stage, a 1 year old recruit stage, and ages 2-10). Lets say there is one stock, with six sub-stocks at unique locations.
### basic spatial properties
n_stages= 10 ### number of stages (eggs, recruits, adults)
n_stocks <- 1 ### number of stocks to simulate
n_loc_stocks <- 6 ### substocks per stock
n_loc = n_stocks*n_loc_stocks ### total number of "stocklets"The maturity of ages 2 and 3 are 25% and 90% followed by 100% thereafter. We can take beverton-holt parameters from the pacific herring stock assessment. We also assume instantaneous natural mortality is 0.58. Assume desnity dependence in recruitment is determined locally, and fish do not follow other fish to spawn (e.g., like MacCall et al. 2018) but this could be done by setting follow to TRUE. Assume no difference in productivity among locations.
maturity = c(0.25,0.9,1,1,1,1,1,1,1) # proportional maturity for ages 2 through plus group
stage_mat <- 2 # stage at first maturity
### taken from Table 13 of the 2015 herring assessment
B0 =26.81/n_loc_stocks # tons in central coast in 2013/14
E0 = B0*1e6*100*1000 # eggs at 6 years old
h_assess <- 0.8 # group steepness
### taken from table E.7 of the 2014 herring assessment
R0 = 287.17e06/n_loc_stocks # maximum local recruitment
beta_0= 1/R0
alpha_0= (E0)/(R0)*(1-h_assess)/(4*h_assess)
dd <- c(1) # density dependence is local
follow <- c(FALSE) # fish do not follow each other (GWOF)
prod_diff = 0.1 # no productivity differential
M_assess <- c(0.58) # natural mortality (fixed)#### stochasticity parmameters
rec_sd <- c(0.8) # log-scale sd in rec-devs
C <- 1e10 ### dirichlet stochasticity in migration (high C is low stochasticity)
obs_sd = c(0.4) ### observation error
h_crit= c(0.25) ### lower biomass threshold for closure
phi <- c(0.6) ### recruitment autocorrelation
surv_rho <-1e-6 ### survival stochasticity (low rho is low stochasticity)### simulation parameters
n_iter = 50 #length of the time series
start = 21 #when to start assessment
warmup = 20 #
warmup2 = 25 # when to start recording data
n_reps= 100 # number of times to run (not run)
s_timeout= 60 # seconds to timeout and restart for a given assessment run
### assessment parameters (bow many samples to take from the simulated population)
ns_samps=200
nf_samps=200### set cauchy scale & recr. corrs to achieve desired retention & syncrony
retention <- function(x)
sapply(x,function(y) mean(diag(spat_cor_mat(n_loc, spat_scale =y, sumto1 = TRUE,offset= 0))))
obj_stray <- 1-c(0.01,0.2,0.4)
obj_juv_stray <- 1-c(0.2,0.4)
ret_fun <- function(x) exp(optim(1,
function(p) (retention(exp(p))-x)^2,
method="BFGS")$par)
stray_s_scale <- sapply(obj_stray,ret_fun)
retention(stray_s_scale)## [1] 0.9900004 0.8000003 0.5999998
juv_stray_s_scale<- sapply(obj_juv_stray,ret_fun)
retention(juv_stray_s_scale)## [1] 0.8000003 0.5999998
offset=0.1
synchrony <- function(n_loc,spat_sd)
mapply(function(y,x) {(1+(y-1)*(mean(spat_cor_mat(y, spat_sd =x)[lower.tri(spat_cor_mat(y, spat_sd =x))])))/y}, n_loc,spat_sd)
obj_synchrony <-c(0.2)
sd_fun <- function(x) exp(optim(2,
function(p) (synchrony(n_loc,exp(p))-x)^2,
method="BFGS")$par)
rec_corr_sd <- sapply(obj_synchrony,sd_fun)
### pre-check for positive definite covariance matrix using cholesky decompose
lapply(rec_corr_sd,function(x) spat_cor_mat(n_loc,spat_sd= x,offset= offset))## [[1]]
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0000000000 0.0987701771 0.0011659042 0.0001266722 0.0011659042
## [2,] 0.0987701771 1.0000000000 0.0987701771 0.0011659042 0.0001266722
## [3,] 0.0011659042 0.0987701771 1.0000000000 0.0987701771 0.0011659042
## [4,] 0.0001266722 0.0011659042 0.0987701771 1.0000000000 0.0987701771
## [5,] 0.0011659042 0.0001266722 0.0011659042 0.0987701771 1.0000000000
## [6,] 0.0987701771 0.0011659042 0.0001266722 0.0011659042 0.0987701771
## [,6]
## [1,] 0.0987701771
## [2,] 0.0011659042
## [3,] 0.0001266722
## [4,] 0.0011659042
## [5,] 0.0987701771
## [6,] 1.0000000000
lapply(rec_corr_sd,function(x) chol(spat_cor_mat(n_loc,spat_sd=x,offset=offset)))## [[1]]
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0.09877018 0.001165904 0.0001266722 1.165904e-03 0.0987701771
## [2,] 0 0.99511027 0.099139787 0.0011590603 1.157224e-05 -0.0086318511
## [3,] 0 0.00000000 0.995072833 0.0991436175 1.169158e-03 0.0008715699
## [4,] 0 0.00000000 0.000000000 0.9950724515 9.914263e-02 0.0010823202
## [5,] 0 0.00000000 0.000000000 0.0000000000 9.950719e-01 0.0990348547
## [6,] 0 0.00000000 0.000000000 0.0000000000 0.000000e+00 0.9901313599
synchrony(n_loc,rec_corr_sd)## [1] 0.1999998
### range of fishery parameters for the simulation
Fmort <- c(0,0.1,0.2,0.25,0.375)
spat_alloc <- c(1)
stage_forecast= c(2)
assessment= TRUENow combine parameter combinations into a dataframe using expand.grid.
### generate the dataframe of parameters and output
param.df <- expand.grid(stray_s_scale= stray_s_scale,
juv_stray_s_scale= juv_stray_s_scale,
rec_corr_sd=rec_corr_sd,
surv_rho= surv_rho,
phi=phi,
Fmort= Fmort,
h_crit= h_crit,
dd= dd,
spat_alloc=spat_alloc,
M_assess=M_assess,
follow= follow,
h= h_assess,
obs_sd= obs_sd,
stage_forecast= stage_forecast,
MPA= c(0,6),
diff= prod_diff)%>%
mutate(a_bh_assess = alpha_0,
b_bh_assess = beta_0,
retention = retention(stray_s_scale),
juv_retention= retention(juv_stray_s_scale),
sync = synchrony(n_loc,rec_corr_sd))
param.df <- subset(param.df,!(follow==TRUE&dd==0)&!(diff!=1&dd==0)&
!(Fmort%in%c(0.025,0.05))&
!(follow==FALSE&juv_stray_s_scale%in%juv_stray_s_scale[c(2,3)])&
!(Fmort==0&h_crit>0.25)&
!(Fmort==0&MPA==6))
param.df[param.df$Fmort==0,"assess"]## NULL
# for sd = 0.7, 6 sites
coefs <- c(1.1180355,-1.0114039,0.3912513,1.4970863,-0.9118604,-0.4836503)
param.df$site_sd <- as.numeric(cbind(1,poly(log10(param.df$rec_corr_sd),5,raw= TRUE))%*%t(t(coefs)))
# set ids by row for later use
param.df$id <- 1:nrow(param.df)We’re using an ADMB assessment, but could avoid an assessment using assessment = FALSE.
### give current directory so we can come back after running admb
wd <- getwd()
#### example simulation function ####
example_sim<- function(x){
setwd(wd)
directory= tempfile(tmpdir= getwd())
dir.create(directory)
file.copy(c("Herring.tpl","Herring","Herring.cpp","Herring.obj","Herring.htp"),directory)
setwd(directory)
rep <- fishery_simulate(n_loc=n_loc, ### total number of subpopulations
n_stages=n_stages, ### number of age/stage classes
stage_mat=2, ### stage/age at maturity
maturity = maturity,
a_bh=1/c(rep(1/x$a_bh_assess*(1-(x$diff-1)/n_loc),5),1/x$a_bh_assess*x$diff), ### BH alpha parameter
b_bh=1/c(rep(1/x$b_bh_assess*(1-(x$diff-1)/n_loc),5),1/x$b_bh_assess*x$diff),
phi=x$phi, ### lag 1 recruitment autocorrelation
M=x$M_assess, ### natural mortality (a constant, vector or matrix)
surv_rho= x$surv_rho, ### random survival stochasticity
assessment= TRUE, ### whether to assess the stock
spat_alloc= 2,
collective_dd=x$dd, ### group density dependence
spat_scale=x$stray_s_scale, ### spatial Cauchy straying scale
juv_spat_scale=x$juv_stray_s_scale, ### spatial Cauchy straying scale
spat_sd=x$rec_corr_sd, ### spatial SD of recruitment synchrony
site_sd=x$site_sd, ### degree of stochasticity in alpha
C=C, ### inverse of degree of stochasticity in straying
obs_sd=x$obs_sd, ### observation error on surveys of total biomass
stock_IDs= rep(1:n_stocks,each= n_loc/n_stocks), ### ids representing to which stock the populations belong
harvest_IDs= c(rep(1,n_loc-1),ifelse(x$MPA==0,1,NA)),
n_iter=n_iter, ### total number of iterations
warmup= warmup,
warmup2= warmup2, ### warmup prior to harvest
point.estimate=FALSE, ### assessment point estimate vs full posterior
h_crit=x$h_crit, ### lower harvest cutoff
Fmort=x$Fmort, ### proportional fishing rate
fit.compile=fit.compile, ### compiled assessment model
ret_ts= TRUE,
follow=x$follow,
stage_forecast= x$stage_forecast)
setwd(wd)
unlink(directory,recursive= TRUE)
return(list(TS= rep$TS,metrics = rep$metrics,rep$age_quantiles))
}data <- rows.to.list(param.df)
out <- example_sim(data[[20]])## iterations 1:20 of 50 complete (warmup)
## iteration 22 of 50 completed 0.257 sec
## iteration 23 of 50 completed 0.224 sec
## iteration 24 of 50 completed 0.242 sec
## iteration 25 of 50 completed 0.234 sec
## iteration 26 of 50 completed 0.264 sec
## iteration 27 of 50 completed 0.273 sec
## iteration 28 of 50 completed 0.283 sec
## iteration 29 of 50 completed 0.3 sec
## iteration 30 of 50 completed 0.306 sec
## iteration 31 of 50 completed 0.284 sec
## iteration 32 of 50 completed 0.312 sec
## iteration 33 of 50 completed 0.294 sec
## iteration 34 of 50 completed 0.262 sec
## iteration 35 of 50 completed 0.3 sec
## iteration 36 of 50 completed 0.316 sec
## iteration 37 of 50 completed 0.346 sec
## iteration 38 of 50 completed 0.346 sec
## iteration 39 of 50 completed 0.377 sec
## iteration 40 of 50 completed 0.366 sec
## iteration 41 of 50 completed 0.398 sec
## iteration 42 of 50 completed 0.432 sec
## iteration 43 of 50 completed 0.417 sec
## iteration 44 of 50 completed 0.511 sec
## iteration 45 of 50 completed 0.506 sec
## iteration 46 of 50 completed 0.448 sec
## iteration 47 of 50 completed 0.47 sec
## iteration 48 of 50 completed 0.5 sec
## iteration 49 of 50 completed 0.521 sec
## iteration 50 of 50 completed 0.537 sec
### format example output
index <-reshape2::melt(out$TS$B_stocks)
names(index) <- c("stock","year","type","index_real")
assess <- melt(out$TS$assess)%>%filter(Var3>warmup&Var4=="Estimated Stock Size")
names(assess) <- c("stock","year","iter","type","est")
assess$type <- "retrospective estimate"
forecast <- melt(out$TS$forecast)%>%filter(!is.na(value)&Var4=="Prediction")
forecast$type <- "forecast (pre-harvest)"
sites <-melt(t(out$TS$B[,,1]),value.name= "biomass")%>%
join(melt(out$TS$Loc_Harvest,value.name= "harvest"))%>%
mutate(Site= factor(Var2,labels= c(1,2,3,4,5,"Closed Area")))
site_mean <- sites%>%
group_by(Var1)%>%
dplyr::summarize(mu= sum(biomass),
muH= sum(harvest))%>%
data.frame()
site_mean$B0_est <- out$TS$B0_est### plot example sims
options <- theme(axis.text.x= element_blank(),
panel.background= element_blank(),
strip.background= element_blank(),
panel.grid.minor= element_blank(),
panel.grid.major.x=element_blank(),
legend.position= c(0.8,0.8),
legend.spacing.y = unit(-0.5, "cm"),
legend.title= element_blank(),
legend.key= element_blank(),
legend.background= element_blank(),
legend.text= element_text(size= 9),
legend.direction= "horizontal",
axis.title= element_text(size= 12),
axis.text.y= element_text(size= 12))
### stock scale plot
stock_plot<- ggplot(aes(year, index_real/out$TS$B0),data=index)+
geom_line(aes(linetype= type,colour= type), size= 1,data= subset(index,type=="actual"))+
geom_line(aes(year,est/out$TS$B0,group= factor(iter),linetype=type,colour= type),
data = subset(assess,type== "retrospective estimate"))+
annotate(x= 10,y= Inf, vjust=1 ,label ="A) stock scale",geom="text")+
geom_segment(aes(y= 0,yend= muH/out$metrics$s_b0
,x= Var1-1,xend= Var1-1),size= 2,col= "grey40",data= site_mean)+
geom_point(aes(Var3,value/out$TS$B0,fill= type),shape= 21,colour= "red",data= forecast)+
geom_point(aes(shape= type),size=1,fill= "grey40",data= subset(index,type== "observed"))+
theme_bw()+
scale_linetype_manual(values= c(1,1),name= "")+
scale_shape_manual(values= c(1),name= "")+
scale_fill_manual(values= c("red","grey40"),name= "")+
scale_colour_manual(values= c("blue",alpha("black",0.4)),name= "")+
geom_line(aes(y= 0.4*B0_est/out$TS$B0,x= Var1-1),linetype= "dotted",col= "red",data= site_mean)+
geom_hline(yintercept= 0.4,col= "red",size= 1,alpha= 0.5)+
ylab("Spawning Biomass")+
#coord_cartesian(ylim= c(0,2),expand= c(0,0),xlim= c(0,n_iter))+
xlab(expression(paste("Time ",symbol('\256'))))+
scale_y_continuous(breaks= c(0,1),labels= c(0,"B0"))+
options+
theme(axis.title.x= element_blank())
### substock scale plot
ages <- t(apply(out$TS$ages,c(1,3),sum))
plot(I(ages[3:50,2]/1e6)~ages[1:48,1],xlim= c(0,max(ages[1:48,1])))
substock_plot <- ggplot(aes(Var1,biomass/out$metrics$p_b0
),data= sites)+
geom_segment(aes(y= 0,yend= muH/out$metrics$s_b0
,x= Var1-1,xend= Var1-1),size= 2,col= "grey40",data= site_mean)+
geom_path(aes(colour= Site,size= Site),lineend= "round")+#facet_grid(Site~.)+
geom_hline(yintercept= 0.25,col= "red",size= 1,alpha= 0.5)+theme_bw()+
ylab("Spawning Biomass")+
annotate(x= 10,y= Inf, vjust=1, label ="B) substock scale",geom="text")+
# coord_cartesian(ylim= c(0,2),expand= c(0,0),xlim= c(0,n_iter))+
xlab(expression(paste("Time ",symbol('\256'))))+
scale_y_continuous(breaks= c(0,1),labels= c(0,"B0"))+
options+
scale_colour_manual(values= c(brewer.pal(7,"PRGn")[c(1:3,6,7)],"black"))+
scale_size_manual(values= c(rep(0.5,5),1))
grid.draw(rbind(ggplotGrob(stock_plot), ggplotGrob(substock_plot), size = "first"))